Anderson Transition in Disordered Bilayer Graphene 
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Employing the Kernel Polynomial method (KPM), we study the electronic properties of the graphene bilayers 
in the presence of diagonal disorder, within the tight-binding approximation. The KPM method enables us to 
calculate local density of states (LDOS) without need to exactly diagonalize the Hamiltonian. We use the 
geometrical averaging of the LDOS's at different lattice sites as a criterion to distinguish the localized states 
from extended ones. We find that bilayer graphene undergoes Anderson metal-insulator transition at a critical 
value of disorder strength. 
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INTRODUCTION 

Graphene (the 2D allotrope of carbon with honeycomb 
structure) has attracted tremendous attraction since its isola- 
tion by Novoselov etal 01- Some of the features that makes 
graphene distinct from previously known materials are proper- 
ties such as, possibility to control the charge carrier types from 
hole to electron in the same sample by a gate voltage, high mo- 
bility of charge carriers (about two orders of magnitude larger 
than the best silicon based semiconductors), anomalous inte- 
ger Quantum Hall effect (IQHE) at room temperature J2, 0] . 
Among them, the high mobility of carriers which is due to 
the vanishing of back-scatterings, makes graphene favorable 
in fabrication of carbon-based electronic devices. 

The spectral and transport properties of graphene are well 
described within the tight-binding approximation. In this 
model the two upper energy bands, which are due to the it 
and 7T* bonds of p z orbitals normal to the honeycomb lattice 
plane, touch each other with a linear dispersion at the corners 
of the Brillouin zone (the so called A' -points) |4j]. In this pic- 
ture graphene is a semi-metal whose low energy excitations 
are massless fermions propagating with Fermi velocity (vf) 
of about 1/300 times the light speed. These excitations are 
called Dirac-fermions, which, along with strong stability of 
graphene structure (which is related to the in plane a bonds), 
are responsible for the odd properties of graphene 10]. 

Employing graphene in electronic devices requires the 
opening of gap in its electronic spectrum. This can be done 
by constraining its geometry either to graphene nano-ribbons 
or quantum-dots. However these methods affect the electronic 
transport because of formation of rough edges and also en- 
hancement of the Coulomb blockade effects in the small size 
structures ^ . Moreover, by reducingthe size and symmetry, 
the selection rules tend to get more restrictive, hence reducing 
the possible conduction channels. Nevertheless, it has been 
shown that by applying an electric field normal to a graphene 
bilayer a gap can be opened whose magnitude is proportional 
to the intensity of the applied electric field la, |7[] . The size of 
gap can be as large as 0.1-0.3 eV, thereby making the bilayer 
graphene an attractive candidate in carbon-based transistors 



The graphene bilayer consists of two layers of carbon atoms 
in honeycomb structure placed on top of each other, making 
the Bernal stacking In this stacking the upper layer is 
rotated 60 degree relative to the lower one, in such a way that 
the atoms in one sublattice, i.e A\ and Ai, are on the top of 
each other while the atoms of the S-sublattice in upper layer 
are placed on the top of the hallow center of hexagons of the 
lower layer . In the absence of external perpendicular elec- 
tric field the band structure of graphene bilayer consists of 
four bands, arising from the 7r-bonding between p z orbitals. 
Two of these bands touch each other at zero energy with a 
quadratic dispersion, and so give rise to the massive quasi- 
particles also called chiral fermions |6J. The other two bands 
are separated by energy scale corresponding to the inter-plane 
hopping energy, t±; one laying below the zero energy and the 
other above it. The electronic structures of bilayer as well 
as multilaye r graph ene have been extensively investigated re- 
cently JlliliMl- 

The single and double layer samples of graphene can be 
grown with high crystalline quality. However there are in- 
evitable sources of disorder which may affect the electronic 
transport in the fabricated samples. Typical forms of disor- 
der includes surface ripples, topological defects, vacancies, 
ad-atoms, charge impurities and polarization field of the sub- 
strate. Scaling theory of localization predicts that all the elec- 
tronic states are localized in two-dimensional systems, once 
smallest amount of disorder is introduced 111 511 . According to 
this theory, the single and bilayer of graphene are expected 
to be insulators at zero temperature. However theoretical as 
well as experimental studies show that both single and bilayer 
graphene, have a minimum conductivity of the order of con- 
ductance quantum (e 2 //i) at the charge neutrality point lfl6l . 
It is also shown that the minimum conductivity is at least twice 
larger in the bilayer graphene lfl4l \l7\\ . These results are in 
contrast to prediction of the scaling theory of localization and 
motivate us to study the localization properties of electronic 
states in graphene systems. In our previous work, we inves- 
tigate the the single layer graphene in the presence of on- 
site disorder by using the numerically powerful approach of 
KPM I1 1811 . There, we found that disordered graphene mono- 
layer shows Anderson transition at a critical value of disorder 
intensity which is of the order of bandwidth. Our calcula- 
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tion also brought an unusual theoretical prediction according 
to which, the localization starts from the states near the Dirac 
points and spreads toward the band edges by increasing the 
amount of disorder 11811 , This type of metal-insulator transi- 
tion driven by short range diagonal disorder due to neutral 
impurities, was experimentally verified in Hydrogen dozed 
graphene by angular resolved photoemission experiment lfl9ll . 

In this paper we investigate the localization properties of 
the electronic states in bilayer graphene in presence of diag- 
onal disorder. We consider minimal coupling between two 
layers of graphene, that is, we take into account only the hop- 
ping between A\ and A2 sublattices. Here we use the kernel 
polynomial method (KPM) lEoll . which consists in the expan- 
sion of various spectral functions in terms of a complete set 
of polynomials. We calculate different local density of states 
(LDOS), from which we introduce a quantity for distinguish- 
ing the localized states from extended ones. The CPU time in 
this method grows proportional to square of the system size, 
enabling us to study large lattice sizes in a moderate time. 



MODEL HAMILTONIAN 

As was mentioned before, we use a minimal tight-binding 
model for describing the low energy transport properties of 
the graphene bilayer. In this model we consider only the hop- 
ping between the p z orbitals residing on the nearest neighbor 
sites. In this picture the two layers are assumed to communi- 
cate only through the hopping among A\ and A2 sublattices 
in the Bernal stacking. The Hamiltonian can be written as: 



W is a measure for the intensity of disorder. Hereafter, we 
assume the unit of energy to be set by t. 



KERNEL POLYNOMIAL METHOD 

Local density of states (LDOS), denoted by pi(E), is a 
quantity that measures the contribution of a given lattice site 
i in total density of energy states of the lattice in the interval 
[E, E + dE], and is defined by the following relation: 
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Oi(E) = J2\{i\E n )\ 2 S(E - E„ 



(2) 
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in which \E n ) is the energy eigenvector corresponding to en- 
ergy eigenvalue E n . \(i\E n )\ 2 in Eq. (f2]i is the probability of 
finding an electron with energy E„ on the site i, which in the 
absence of disorder is the same for all lattice sites as a result of 
translational invariance. However for localized states encoun- 
tered in disordered systems, this probability drastically varies 
over the system. If the Hamiltonian possesses an extended 
eigenstate with the eigenvalue between E and E + dE, then 
all sites are comparably likely to be present in this state, while 
for a localized state only a limited number of sites have appre- 
ciable probability of being occupied. Therefore LDOS would 
be a suitable quantity to distinguish an extended state from a 
localized one. This can be accomplished by comparing the 
geometric and arithmetic averagings of LDOS's at different 
lattice points. The geometric average of LDOS's known as 
typical DOS is defined as: 



H= - tJ2J2 a l,i b m,j + h.c. (1) 
- t± a\ A a 2 ,i + h.c. 

i 
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+ ^ ] ^ t e m,i a m,i a "m,i + e rn.i^m,i^"i,i, 
m—1 i 

in which a' m j(a TO ,j) creates (annihilates) a ^-electron at site 
(m,i) on the sublattice A, with m = 1,2, where m is the 
index of layers and i labels the sites on each A-sublattice in 
given layer. Similarly, b* m i and b m _i are the corresponding 
creation and annihilation operators in the B-sublattices of the 
two layers. Here, t denotes the intra-plane hopping integral 
between the nearest neighbors and t± is the inter-plane hop- 
ping amplitude between two layers. Empirical estimates of the 
hopping terms are t ~ 3.16 eV, and t± ~ 0.39 eV 112 ll 12211 . 



There are further inter-plane couplings, where for simplicity 
we ignore them in our study. Also, i and e h m i in the last 
term denote the on-site energies at A and B sublattices on each 
layer, respectively. To introduce disorder on the model, we 
choose the on-site energies randomly from a uniform distribu- 
tion in the interval [— This is the so called Anderson 
model with spatially uncorrected diagonal disorder in which 



Ptyp{E) = exp 
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K r K, 



(3) 



where K s is the number of sites in a given realization that 
LDOS is calculated and K r is the number of realizations. On 
the other hand, the total density of states is obtained by sum- 
ming up the partial LDOS of all lattice sites, which amounts 
to the following arithmetic average: 



K r K 



D-l 



Pav(E) 



K r K 



where D is the dimension of Hilbert space of the Hamiltonian. 
For an extended state p typ (E) as p av (E), while in the case of 
localized states p typ (E) <C p av (E). 

LDOS is a site dependent quantity which we use KPM [20] 
to compute it. The basic idea of KPM is to expand the spectral 
functions such as, Pi(E), in terms of orthogonal polynomials 
of energy E. In principle, one can use any kind of orthogonal 
polynomials. In this paper we use Chebyshev polynomials. 
Therefore we expand LDOS as: 



Pi(E) 



ttVI - E 2 
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Po 



2j2PnT n {E) 



(5) 
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where the coefficient /x„ is given by 12011 : 



Pi{E)T n (E)dE = ^(i\T n (H)\i). (6) 



in the above equation is the rescaled Hamiltonian which 
can be obtained by a simple shift and scaling transformation 
of the original Hamiltonian to ensure that eigenvalues of H lay 
in the interval [—1,1]. A similar procedure for the expansion 
of total DOS gives the following moments: 



Pav{E)T n (E)dE = -Tr[T n (H) 



(7) 



The moments given in Eqs. (O and ©, can be evaluated em- 
ploying the recursion relation of Chebyshev polynomials first 
discussed by Wang |23|]. The matrix elements are computed 
on the fly, without saving any matrix, which is a key aspect of 
KPM. The summation over the diagonal matrix elements re- 
quired in the trace in Eq. ^} can be performed stochastically 
which dramatically reduces the computer time. These points 
makes it possible to study reasonably large systems even on a 
desktop computer lEoll . 

When the series expansion of Eq. (O is used in numeri- 
cal calculations, it has to be truncated at a finite order N. 
This truncation leads to the infamous Gibbs oscillations in the 
LDOS's. To relieve this effect, some standard damp ing fac- 
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tors 112411 (called g-factors) have been suggested [2 
using which Eq. (O can be evaluated as follows: 



Pi{E) 



1 



Wl - E 2 



N 
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(8) 



In our work, we found that the following positive g-factor 
(Jackson's) 

(iV-n + l)cos( 1 ^ T ) + S in( 1 ^ T )cot( 1 ^ T ) 

9n = , (9) 

is suitable to damp the oscillations arose in the calculation of 
LDOS's H. 



RESULTS AND DISCUSSIONS 

We operate with the Hamiltonian given in Eq. (Q~|i with dif- 
ferent strengths of disorder (W) to investigate the Anderson 
transition on the bilayer graphene and then calculate LDOS's 
by means of KPM. We carried out the calculations for the lat- 
tices consisting of L = 20 x 10 3 to L = 80 x 10 3 sites. In order 
to obtain reliable results, the order of expansion, N, has to be 
reapportioned depending on the system size as N ~ L 2 . The 
number of random lattice sites used in averaging the LDOS's 
is K s = 15 for each of the K r = 15 different realizations. 
In Fig. [D(a) we compare the total and typical DOS corre- 
sponding to L = 8 x 10 3 sites and N = 4000 moments for 
W/t = 0.5 and W/t = 1.0. As can be seen in this figure 



Pty P is nonzero and almost equal to p av for the entire energy 
bandwidth, indicating that none of the states are localized for 
these strengths of disorder. The four sharp peaks located near 
E/t~ ±1.0 in the W/t = 0.5 plot are the Van-Hove singular- 
ities due to the four saddle points in the band structure of the 
clean bilayer graphene. Therefore, this disorder strength does 
not appreciably alter the overall spectral aspects of the clean 
bilayer graphene. There are also four jumps near E/t = ±3, 
corresponding to the four extrema in energy surfaces. As can 
be seen in Fig.[T]-(b), for W/t > 1.0 these singularities tend 
to smear out. One has to note that the total bandwidth of clean 
bilayer graphene is very close to the value 3t of a mono-layer 
sample, due to the dominant in-plane energy scale, t. How- 
ever, as we will see shortly, when the disorder is introduced, 
the presence of the second layer causes the states in the bi- 
layer graphene resist more against Anderson localization than 
the mono-layer. 

The results for the larger values of the disorder strength are 
displayed in Fig. [2] As can be seen in this figure, the aver- 
age density of states still resembles that of the clean graphene 
bilayer up to W/t = 2. Beyond W/t ~ 2, the spectral prop- 
erties start to significantly deviate from the clean bilayer. The 
typical DOS starts to vanish at two bounds of the energy spec- 
trum. This behavior is quite similar to the localization be- 
havior of 3D bands, where a mobility edge sets in, rendering 
the states at the tails localized. Further increase in the disor- 
der strength, eventually localizes the entire spectrum beyond 
W / 1 ~ 8. The fact that states at the band edge are more sensi- 
tive to disorder and get localized before those at the center of 
the band is in bold contrast to mono-layer graphene IU8U25 I. 
This unusual behavior of mono-layer graphene has been ex- 
perimentally observed in photoemission experiment lfl9ll . 

Fig.[3]shows intensity plot of LDOS in the xy-plane of the 
lattice for different values of disorder near charge neutrality 
point, E = 0. As can be seen, by increasing the disorder 
strength, the spatial distribution of the DOS tends to clump 
into clusters. At the disorder strength W/t = 8.0, these clus- 
ters of non-zero LDOS become entirely disconnected in such 
a way that the energy states around E = would no longer 
contribute to the electric conduction. 

To rule out the possibility of finite size artifacts, we perform 
a finite size scaling analysis. Let us define R(E) as the ratio 
of the typical to average DOS for a given electronic mode at 
energy E 02611 : 



R(E) = 



PtypjE) 
Pav(E) 



(10) 



Since always the arithmetic averages of positive numbers are 
greater than the geometric average, one always has, 



R(E) < 1, ME. 



(11) 



In the absence of disorder the equality is realized, while turn- 
ing the disorder on, reduces R(E) from 1. We computed this 
quantity for different values of disorder strength and different 
lattice sizes. The results for two energies E = and E = 0.2 
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are depicted in Figs. [4] and [5] respectively. The lattice sizes 
are L = 2 x 10 4 , 4 x 10 4 , 6 x 10 4 and 8 x 10 4 . As can be 
seen in a wide range of disorder strengths, this ratio tends to 
converge almost to the same curve with increasing L, which 
indicates the validity of our results on the infinite lattice size. 
We have checked that increasing the order of expansion from 
N = 4000 to TV = 12000, does not appreciably change our 
results. 

CONCLUSION 

In summary, using KPM method to find the local density 
of states; and studying their geometrical averages, we inves- 
tigated the effect of uncorrelated disorder on the electronic 
states of the graphene bilayer. The localization behavior of 
bilayer graphene is reminiscent of 3D bands, i.e. the localiza- 
tion starts from the states at the edges of the energy spectrum. 
The states near the charge neutrality point E = remains ex- 
tended up to very large value of disorder strength W/ 1 ~ 8. 
Therefore, our finding shows that the bilayer graphene re- 
mains metal in the presence of on-site disorder up to a crit- 
ical value of disorder at which undergoes the Anderson metal- 
insulator transition. This critical value is about three times 
larger than the critical value that we obtained for Anderson 
transition in single-layer graphene flUl . The inter-layer hop- 
ping seems to be responsible for such a drastic difference be- 
tween the localization properties of mono-layer and bilayer 
graphene. This difference manifests itself in more than twice 
a large minimal conductivity of bilayer graphene compared to 
mono-layer samples. The precise understanding of how the 
(small) inter-layer hopping can lead to such a remarkable dif- 
ference in the localization properties remains an open question 
which needs further investigations. 
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FIG. 1 : (Color online) Average (red) and typical (black) DOS for the 
lattice size L = 80 x 10 3 with K r x K s = 15 x 15 and N = 4000 
at disorder intensities (a) W/t = 0.5 and (b) W/t = 1.0 
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FIG. 2: (Color online) Average (red) and typical (black) DOS for the 
lattice size L = 80 x 10 3 with if r x K s = 15 x 15 and N = 4000 at 
disorder intensities (a) W/t = 2.0, (b) W/t = 4.0, (c) W/t = 6.0 
and (d) W/t = 8.0. 
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FIG. 3: (Color online) LDOS map for various disorder intensities for 
E 6 [0.0, 0.01]. In the weak disorder regime, the electron density 
is almost uniformly distributed over the entire lattice. By increasing 
W/t to 8, the electron density becomes confined in disconnected 
regions of lattice. 
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FIG. 4: (Color online) The ratio of typical to average DOS (R(E)), 
computed at the band center (E — 0), versus the strength of disorder 
for different lattice sizes L = 2 x 10 4 , 4 x 10 4 , 6 x 10 4 , 8 x 10 4 . 
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FIG. 5: (Color online) The ratio of typical to average DOS (R(E)), 
computed at E — 0.2 , versus the strength of disorder for different 
lattice sizes L = 2 x 10 4 , 4 x 10 4 , 6 x 10 4 , 8 x 10 4 . 



